{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "arabic-friendly",
   "metadata": {},
   "source": [
    "# 频率分析 (3)：热力学能矫正"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exciting-twelve",
   "metadata": {},
   "source": [
    "> 创建时间：2021-06-18"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fatty-shoulder",
   "metadata": {},
   "source": [
    "在这份文档中，我们将简单地讨论从 Gaussian 生成的 formated checkpoint 文件 (fchk 或 fch 后缀名)，产生热力学能量矫正。\n",
    "\n",
    "我们计算的分子与 [频率分析 (1)](freq_1) 相同，为没有优化到能量最低结构的 C<sub>2</sub>O<sub>4</sub>H<sup>+</sup> 离子。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interim-arrow",
   "metadata": {},
   "source": [
    "该文档的重要的参考资料是 Gaussian 的白皮书 Thermochemistry in Gaussian [^gaussian-thermo]。我们所使用的公式也 (**很无耻地**) 会与该文档几乎完全相同；并且作为程序实现笔记，也不会对具体的公式推导作讨论。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "immune-railway",
   "metadata": {},
   "source": [
    ":::{warning}\n",
    "\n",
    "不处于能量最低结构的分子一般来说不适合用作频率分析。此时不仅分子光谱、同时热力学矫正从理论上也是不允许的。\n",
    "\n",
    "这份文档尽管使用了有虚频的分子，但若要进行有价值的热力学矫正，仍然需要先对分子的结构进行优化。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mighty-neighborhood",
   "metadata": {},
   "source": [
    "分子对应的输入卡 {download}`C2O4H.gjf`、输出文件 {download}`C2O4H.out` 与 fchk 文件 {download}`C2O4H.fchk` 在链接中可供下载。这份文档的目标将是重复输出文件中的热力学矫正量。其一部分输出是："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "broken-bulgaria",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Thermal correction to Energy=                    0.030349\n",
      " Thermal correction to Enthalpy=                  0.031293\n",
      " Thermal correction to Gibbs Free Energy=        -0.002096\n",
      " Sum of electronic and zero-point Energies=           -370.093195\n",
      " Sum of electronic and thermal Energies=              -370.088816\n",
      " Sum of electronic and thermal Enthalpies=            -370.087872\n",
      " Sum of electronic and thermal Free Energies=         -370.121261\n",
      " \n",
      "                     E (Thermal)             CV                S\n",
      "                      KCal/Mol        Cal/Mol-Kelvin    Cal/Mol-Kelvin\n",
      " Total                   19.044             14.495             70.273\n",
      " Electronic               0.000              0.000              0.000\n",
      " Translational            0.889              2.981             39.370\n",
      " Rotational               0.889              2.981             26.171\n",
      " Vibrational             17.267              8.534              4.732\n",
      " Vibration     1          0.683              1.701              1.500\n",
      " Vibration     2          0.731              1.565              1.145\n",
      " Vibration     3          0.897              1.158              0.562\n",
      " Vibration     4          0.941              1.067              0.479\n",
      "                       Q            Log10(Q)             Ln(Q)\n",
      " Total Bot       0.920866D+01          0.964196          2.220144\n",
      " Total V=0       0.811762D+13         12.909429         29.725059\n",
      " Vib (Bot)       0.238627D-11        -11.622280        -26.761289\n"
     ]
    }
   ],
   "source": [
    "with open(\"C2O4H.out\", \"r\") as f:\n",
    "    while \"Zero-point correction\" not in f.readline(): continue\n",
    "    for _ in range(23): print(f.readline()[:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recreational-sister",
   "metadata": {},
   "source": [
    "我们的文档的内容原则上可以重现所有的上述数据。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "polish-examination",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "热力学能量矫正本质上是统计热力学的初步应用。我们这里遇到的分子是最为常见的没有对称性、单重态的分子。但一些特殊的情况，譬如具有对称性、多重态等情况，需要对下述代码作改动。这些改动我们不在文档中作补充，因此读者还是需要参考 Gaussian 白皮书 [^gaussian-thermo]。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stopped-marine",
   "metadata": {},
   "source": [
    "## 准备工作"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hired-camel",
   "metadata": {},
   "source": [
    "频率分析部分是由 {download}`freqanal.py` 完成的；它还需要调用读取 Gaussian formchk 文件的小程序 {download}`formchk_interface.py`。这些程序可以下载。频率分析的具体做法已经在 [频率分析 (1)](freq_1) 有所陈述；对于比较一般的非线性分子，它应当可以输出与 Gaussian 近乎一致的频率结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "portable-gibson",
   "metadata": {},
   "outputs": [],
   "source": [
    "from freqanal import FreqAnal\n",
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(5, linewidth=150, suppress=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "labeled-functionality",
   "metadata": {},
   "source": [
    "为了能进行单位换算，我们还需要定义一些常数。这里对它们在这份文档中的变量名与符号作说明，并给出大致的数量。\n",
    "\n",
    "- `E_h` $E_\\mathrm{H}$：Hartree 能量 $4.360 \\times 10^{-18} \\ \\mathrm{J}$；\n",
    "- `a_0` $a_0$：Bohr 半径 $5.292 \\times 10^{-11} \\ \\mathrm{m}$\n",
    "- `N_A` $N_\\mathrm{A}$：Avogadro 常数 $6.022 \\times 10^{23} \\ \\mathrm{mol}^{-1}$\n",
    "- `c_0` $c$：真空光速 $2.998 \\times 10^{8} \\ \\mathrm{m} \\ \\mathrm{s}^{-1}$\n",
    "- `k_B` $k_\\mathrm{B}$：Boltzmann 常数 $1.381 \\times 10^{-23} \\ \\mathrm{J} \\ \\mathrm{K}^{-1}$\n",
    "- `R` $R$：Mole 气体常数 $8.314 \\ \\mathrm{J} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$\n",
    "- `h` $h$：Planck 常数 $6.626 \\times 10^{-34} \\ \\times{J} \\ \\mathrm{s}$\n",
    "- `P_0` $P_0$：标准大气压 $101325. \\mathrm{kg} \\ \\mathrm{m}^{-1} \\mathrm{s}^{-2}$\n",
    "- `amu` $m_\\mathrm{u}$：原子质量单位 $1.661 \\times 10^{-27} \\ \\mathrm{kg}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "productive-activity",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://docs.scipy.org/doc/scipy/reference/constants.html\n",
    "from scipy import constants\n",
    "from scipy.constants import physical_constants\n",
    "\n",
    "E_h = physical_constants[\"Hartree energy\"][0]\n",
    "a_0 = physical_constants[\"Bohr radius\"][0]\n",
    "N_A = physical_constants[\"Avogadro constant\"][0]\n",
    "c_0 = physical_constants[\"speed of light in vacuum\"][0]\n",
    "k_B = physical_constants[\"Boltzmann constant\"][0]\n",
    "R = physical_constants[\"molar gas constant\"][0]\n",
    "h = physical_constants[\"Planck constant\"][0]\n",
    "P_0 = physical_constants[\"standard atmosphere\"][0]\n",
    "amu = physical_constants[\"atomic mass constant\"][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "black-debut",
   "metadata": {},
   "source": [
    "- Calorie 与 Joule 的换算比例 4.184\n",
    "- `pi` $\\pi$：圆周率"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "premium-debate",
   "metadata": {},
   "outputs": [],
   "source": [
    "calorie = 4.184\n",
    "pi = np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gross-delay",
   "metadata": {},
   "source": [
    "我们在整个文档中，使用的温度是 $T = 298.15 \\ \\mathrm{K}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "incorporated-shuttle",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 298.15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pediatric-snake",
   "metadata": {},
   "source": [
    "对 C<sub>2</sub>O<sub>4</sub>H<sup>+</sup> 离子的频率分析对象储存在变量 `fa` 中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "demonstrated-difference",
   "metadata": {},
   "outputs": [],
   "source": [
    "fa = FreqAnal(\"C2O4H.fchk\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "suburban-veteran",
   "metadata": {},
   "source": [
    "整个热力学分析过程需要分为平动 (translation)、电子态 (electronic)、转动 (rotation)、振动 (vibrational) 四部分。需要计算的基础热力学矫正量是熵 $S$ (entropy)、内能 $E$ (thermo energy)、热容 $C$ (heat capacity)。这里的热容是指恒容热容。在计算熵时，我们还需要给出配分函数 $q$ (partition function)。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "configured-retro",
   "metadata": {},
   "source": [
    "## 平动 (Translation)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exciting-treasury",
   "metadata": {},
   "source": [
    "- 配分函数\n",
    "\n",
    "    $$\n",
    "    q_\\mathrm{t} = \\left( \\frac{2 \\pi m k_\\mathrm{B} T}{h^2} \\right)^{3/2} \\frac{k_\\mathrm{B} T}{P_0}\n",
    "    $$\n",
    "    \n",
    "    `m` $m$ 指分子质量，Gaussian 的输出是原子质量单位。它需要先转为 SI 单位制再进行计算。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "stock-president",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = fa.mol_weight.sum() * amu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "institutional-artist",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "32994942.55727539"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_t = (2 * pi * m * k_B * T / h**2)**(3/2) * (k_B * T / P_0)\n",
    "q_t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "premier-concern",
   "metadata": {},
   "source": [
    "- 熵 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    S_\\mathrm{t} = R (\\ln q_\\mathrm{t} + 1 + 3/2)\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "light-encoding",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39.37022220447884"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_t = R * (np.log(q_t) + 1 + 3/2) / calorie\n",
    "S_t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "revolutionary-knife",
   "metadata": {},
   "source": [
    "- 内能 ($\\mathrm{kcal} \\ \\mathrm{mol}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    E_\\mathrm{t} = \\frac{3}{2} R T\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "satellite-student",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8887274245542662"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_t = 3/2 * R * T / 1000 / calorie\n",
    "E_t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "decreased-contribution",
   "metadata": {},
   "source": [
    "- 热容 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    C_\\mathrm{t} = \\frac{3}{2} R\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "executed-section",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.9808063879063096"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_t = 3/2 * R / calorie\n",
    "C_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "neural-latitude",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "functional-leather",
   "metadata": {},
   "source": [
    "## 电子态 (Electronic)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "respected-traveler",
   "metadata": {},
   "source": [
    ":::{warning}\n",
    "\n",
    "电子态的配分函数 $q_\\mathrm{e}$ 受分子多重度的影响；这也会同时影响到熵矫正 $S_\\mathrm{e}$。因此对于多重度不为 1 的分子，下述代码将需要作修改。具体地来说，需要将 `ω_0` $\\omega_0$ 改为多重度的数值。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "strategic-thirty",
   "metadata": {},
   "source": [
    "- 配分函数\n",
    "\n",
    "    $$\n",
    "    q_\\mathrm{e} = \\omega_0\n",
    "    $$\n",
    "    \n",
    "    `ω_0` $\\omega_0$ 指分子多重度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "communist-numbers",
   "metadata": {},
   "outputs": [],
   "source": [
    "ω_0 = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "above-specialist",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_e = ω_0\n",
    "q_e"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "developmental-approach",
   "metadata": {},
   "source": [
    "- 熵 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    S_\\mathrm{e} = R \\ln q_\\mathrm{e}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "personalized-forward",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_e = R * np.log(q_e) / calorie\n",
    "S_e"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "capital-airline",
   "metadata": {},
   "source": [
    "- 内能 $E_\\mathrm{e}$ 与热容 $C_\\mathrm{e}$ 取零值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "final-formula",
   "metadata": {},
   "outputs": [],
   "source": [
    "E_e = C_e = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cubic-corpus",
   "metadata": {},
   "source": [
    "## 转动 (Rotation)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "usual-class",
   "metadata": {},
   "source": [
    ":::{warning}\n",
    "\n",
    "分子转动的配分函数、熵、内能与热容计算都受分子本身构型而影响。\n",
    "\n",
    "- 我们这里计算的是不具有对称性的分子；\n",
    "\n",
    "- 对于线性分子，所有物理量的计算都将发生变化，**这里的代码将不能使用**；\n",
    "\n",
    "- 对于具有对称性的分子需要修改代码；具体地来说，需要将 `σ_r` $\\sigma_\\mathrm{r}$ 设置为分子的对称数；对称数应当指分子的一些原子置换，但通过一系列旋转或反映操作可以重现原来分子的总数量，对于水为 2、氨气为 3、甲烷或苯为 6。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "quality-clinic",
   "metadata": {},
   "outputs": [],
   "source": [
    "q_r = 1.16957e5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "respiratory-participation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "26.17060894487201"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_r = R * (np.log(q_r) + 3/2) / calorie\n",
    "S_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "endangered-pitch",
   "metadata": {},
   "source": [
    "- 转动惯量 `Ixyz` $I_x, I_y, I_z$ ($\\mathrm{kg} \\ \\mathrm{m}^2$)。需要注意，这里并非真的是绕 $x, y, z$ 轴转动的惯量，即\n",
    "    \n",
    "    $$\n",
    "    I_x \\neq \\sum_{A} m_\\mathrm{A} r_{Ax}^2 \n",
    "    $$\n",
    "    \n",
    "    它需要通过对 $3 \\times 3$ 的转动惯量矩阵 $I_{xx}, I_{xy}, \\cdots, I_{zz}$ 作对角化得到。对角化同时会得到三个转动主轴；由于这三个转动主轴相互垂直，确实地可以构建坐标系，因此才称为 $I_x, I_y, I_z$，但这里的 $x, y, z$ 相对于输入卡的坐标系一般有旋转。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "instant-pregnancy",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.33218e-45, 2.34564e-45, 3.43473e-45])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ixyz = fa.rot_eig * amu * a_0**2\n",
    "Ixyz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "related-worship",
   "metadata": {},
   "source": [
    "- 转动特性温度 `Θxyz_r` $\\Theta_{\\mathrm{r}, x}, \\Theta_{\\mathrm{r}, y}, \\Theta_{\\mathrm{r}, z}$ ($\\mathrm{K}$)\n",
    "\n",
    "    $$\n",
    "    \\Theta_{r, x} = \\frac{h^2}{8 \\pi^2 I_x k_\\mathrm{B}}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "collaborative-brunswick",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.30233, 0.1717 , 0.11726])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Θxyz_r = h**2 / (8 * pi**2 * Ixyz * k_B)\n",
    "Θxyz_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "studied-quick",
   "metadata": {},
   "source": [
    "- 配分函数\n",
    "\n",
    "    $$\n",
    "    q_\\mathrm{r} = \\frac{\\pi^{1/2}}{\\sigma_r} \\left( \\frac{T^3}{\\Theta_{\\mathrm{r}, x} \\Theta_{\\mathrm{r}, y} \\Theta_{\\mathrm{r}, z}} \\right)^{1/2}\n",
    "    $$\n",
    "    \n",
    "    其中 $\\sigma_\\mathrm{r}$ 为分子的对称数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "broad-china",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "116957.22677656508"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "σ_r = 1\n",
    "q_r = (pi**(1/2) / σ_r) * (T**3 / Θxyz_r.prod())**(1/2)\n",
    "q_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rubber-limit",
   "metadata": {},
   "source": [
    "- 熵 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    S_\\mathrm{r} = R (\\ln q_\\mathrm{r} + 3/2)\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "stunning-surfing",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "26.170612798005376"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_r = R * (np.log(q_r) + 3/2) / calorie\n",
    "S_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olympic-martin",
   "metadata": {},
   "source": [
    "- 内能 ($\\mathrm{kcal} \\ \\mathrm{mol}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    E_\\mathrm{r} = \\frac{3}{2} R T\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "beautiful-helena",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8887274245542662"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_r = 3/2 * R * 298.15 / 1000 / calorie\n",
    "E_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "virtual-twist",
   "metadata": {},
   "source": [
    "- 热容 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    C_\\mathrm{r} = \\frac{3}{2} R\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "fewer-confidence",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.9808063879063096"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_r = 3/2 * R / calorie\n",
    "C_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dated-virtue",
   "metadata": {},
   "source": [
    "## 振动 (Vibrational)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "analyzed-horizon",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "热力学问题的一般都要求假设分子处于稳定构型 (因此才能有热力学稳定的状态，各种热力学能量作为温度的状态函数才能成立)。因此，具有振动虚频的分子一般认为是**不能**进行振动分析的。\n",
    "\n",
    "我们这里的计算单纯地是 Gaussian 的结果作比对。Gaussian 在计算热力学能时，对虚频的处理是直接忽视。\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "harmful-notification",
   "metadata": {},
   "source": [
    "- 振动特征温度 ($\\mathrm{K}$)\n",
    "\n",
    "    $$\n",
    "    \\Theta_{\\mathrm{v}, K} = h \\nu_K = \\varpi_K h c\n",
    "    $$\n",
    "    \n",
    "    其中，$K$ 是指震动模式，$\\varpi_K$ 是以长度为量纲的振动频率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "plain-communication",
   "metadata": {},
   "outputs": [],
   "source": [
    "ΘK_v = fa.freq[fa.freq > 0] / 1e-2 * h * c_0 / k_B"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "packed-therapy",
   "metadata": {},
   "source": [
    "- 配分函数\n",
    "\n",
    "    $$\n",
    "    q_\\mathrm{v} = \\prod_K \\frac{e^{- \\Theta_{\\mathrm{v}, K} / 2 T}}{1 - e^{- \\Theta_{\\mathrm{v}, K} / T}}\n",
    "    $$\n",
    "    \n",
    "    这里采用的是与 Gaussian 的最终热力学矫正一致的输出 (即 `BOT` 结果)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "incorporated-quebec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.386188575409033e-12"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qK_v = np.exp(-ΘK_v / (2 * T)) / (1 - np.exp(-ΘK_v / T))\n",
    "q_v = qK_v.prod()\n",
    "q_v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recovered-johnston",
   "metadata": {},
   "source": [
    "- 熵 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    S_\\mathrm{v} = R \\sum_K \\left( \\frac{\\Theta_{\\mathrm{v}, K} / T}{e^{\\Theta_{\\mathrm{v}, K} / T} - 1} - \\ln (1 - e^{- \\Theta_{\\mathrm{v}, K} / T}) \\right)\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "narrative-stretch",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.732478685230292"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SK_v = R * (ΘK_v / T / (np.exp(ΘK_v / T) - 1) - np.log(1 - np.exp(-ΘK_v / T))) / calorie\n",
    "S_v = SK_v.sum()\n",
    "S_v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "clear-scheme",
   "metadata": {},
   "source": [
    "- 内能 ($\\mathrm{kcal} \\ \\mathrm{mol}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    E_\\mathrm{v} = R \\sum_K \\Theta_{\\mathrm{v}, K} \\left( \\frac{1}{2} + \\frac{1}{e^{\\Theta_{\\mathrm{v}, K} / T} - 1} \\right)\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "quiet-caution",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17.266670082670938"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EK_v = R * (ΘK_v * (1/2 + 1 / (np.exp(ΘK_v / T) - 1))) / 1000 / calorie\n",
    "E_v = EK_v.sum()\n",
    "E_v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alien-archive",
   "metadata": {},
   "source": [
    "- 热容 ($\\mathrm{cal} \\ \\mathrm{mol}^{-1} \\ \\mathrm{K}^{-1}$)\n",
    "\n",
    "    $$\n",
    "    C_\\mathrm{v} = R \\sum_K e^{- \\Theta_{\\mathrm{v}, K} / T} \\left( \\frac{\\Theta_{\\mathrm{v}, K} / T}{e^{- \\Theta_{\\mathrm{v}, K} / T} - 1} \\right)^2\n",
    "    $$\n",
    "    \n",
    "    这里 Gaussian 白皮书的记号有略微的错误。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "lined-politics",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.533482711830887"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CK_v = R * (np.exp(-ΘK_v / T) * (ΘK_v / T / (np.exp(-ΘK_v / T) - 1))**2) / calorie\n",
    "C_v = CK_v.sum()\n",
    "C_v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stylish-wagon",
   "metadata": {},
   "source": [
    "## 总热力学矫正量"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hidden-category",
   "metadata": {},
   "source": [
    "最终的热力学矫正量直接通过将四部分 (平动、电子态、转动、振动) 相加即可。对于配分函数，其结果通过相乘得到。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "alpine-huntington",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.04412493177947"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_thermal = E_t + E_e + E_r + E_v\n",
    "E_thermal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "equipped-holder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.495095487643507"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_thermal = C_t + C_e + C_r + C_v\n",
    "C_thermal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "configured-order",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70.2733136877145"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_thermal = S_t + S_e + S_r + S_v\n",
    "S_thermal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "separate-monaco",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.208294504188077"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_total = q_t * q_e * q_r * q_v\n",
    "q_total"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "printable-cardiff",
   "metadata": {},
   "source": [
    "## 最终热力学矫正"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worst-attribute",
   "metadata": {},
   "source": [
    "化学中关心的通常是零点能 (zero point energy, ZPE)、焓 (Enthalpy)、Gibbs 自由能 (Gibbs free energy)。这些量不难通过上面的结果给出。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "humanitarian-session",
   "metadata": {},
   "source": [
    "- 零点能 (zero point energy, ZPE) ($E_\\mathrm{H}$, Hartree)\n",
    "\n",
    "    $$\n",
    "    E_\\mathrm{ZPE} = \\sum_{K} \\frac{1}{2} \\Theta_{\\mathrm{v}, K} k_\\mathrm{B}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "hourly-welding",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.025969755181835332"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_zpe = ΘK_v.sum() / 2 * k_B / E_h\n",
    "corr_zpe"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olympic-london",
   "metadata": {},
   "source": [
    "- 内能矫正 ($E_\\mathrm{H}$, Hartree) 实质上就是刚才结果的单位变换而已。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "liable-asbestos",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03034874486989149"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_thermal = E_thermal / (E_h * N_A / 1000 / calorie)\n",
    "corr_thermal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "laughing-percentage",
   "metadata": {},
   "source": [
    "- 焓 (Enthalpy) 矫正 ($E_\\mathrm{H}$, Hartree)\n",
    "\n",
    "    $$\n",
    "    H_\\mathrm{corr} = E_\\mathrm{corr} + k_\\mathrm{B} T\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "prepared-buffer",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03129292973753578"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_enthalpy = corr_thermal + k_B * T / E_h\n",
    "corr_enthalpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "collected-assault",
   "metadata": {},
   "source": [
    "- Gibbs 自由能 (Gibbs free energy) 矫正 ($E_\\mathrm{H}$, Hartree)\n",
    "\n",
    "    $$\n",
    "    G_\\mathrm{corr} = H_\\mathrm{corr} - T S_\\mathrm{tot}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "nonprofit-entity",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.002096189219235073"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_gibbs = corr_enthalpy - T * S_thermal / (E_h * N_A / calorie)\n",
    "corr_gibbs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "anticipated-accused",
   "metadata": {},
   "source": [
    "[^gaussian-thermo]: <http://gaussian.com/thermo/>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
